function [x,fval,residual,exitflag,output]=optim_fitsurface_tau(Param,xdata,ydata,yerr,LB,UB)
options=optimset('MaxIter',500,'TolFun',1e-8,'Display','iter','TolX',1e-8,'MaxFunEvals',2000);
[x,fval,residual,exitflag,output]=lsqcurvefit(@nestedfun,Param,xdata,(ydata./yerr),LB,UB,options);
    function y=nestedfun(Param,xdata)
        C=Param(1);
        qc=Param(2);
        Z=Param(3);
        tau=Param(4);
        qtau=Param(5);
        
        qz=xdata;
        Energy=13.474;
        lambda=0.9201;
        Attl=4854.22e4; % Attenuation length in unit of Angstrom for water @ 13.474 keV
        mu=1/Attl;
        beta=mu*lambda/4/pi;
        
        D=fitsubphase(qz,qc,Energy,Attl);
        y=C*Transmission(qz,qc,beta,lambda).*exp(-Z./D).*exp(-0.5*(1+erf((qz-qtau)/0.0001))*tau.*(qz-qtau)./qtau);
        y=(y./yerr);
    end
end